Numerical study of three-body recombination 
for systems with many bound states 
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Three-body recombination processes are treated numerically for a system of three identical bosons. 
The two-body model potentials utilized support many bound states, a major leap in complexity that 
produces an intricate structure of sharp nonadiabatic avoided crossings in the three-body hyperra- 
dial adiabatic potentials at short distances. This model thus displays the usual difficulties of more 
realistic systems. To overcome the numerical challenges associated with sharp avoided crossings, 
the slow variable discretization (SVD) approach is adopted in the region of small hyperradii. At 
larger hyperradii, where the adiabatic potentials and couplings are smooth, the traditional adia- 
batic method suffices. Despite the high degree of complexity, recombination into deeply bound 
states behaves regularly due to the dominance of one decay pathway, resulting from the strong 
coupling between different recombination channels. Moreover, the usual Wigner threshold law must 
be modified for excited incident recombination channels. 



I. INTRODUCTION 



Three-body recombination has attracted much theo- 
retical and experimental research interest in recent years. 
Recombination is the process in which three free particles 
collide to form a two-body state and a free particle, with 
the released kinetic energy being distributed between the 
final collisional partners. Such reactions are common and 
important in chemical reactions and in atomic, molecu- 
lar, and nuclear physics. In ultracold degenerate Fermi 
gases [l| recombination has been used as a process to 
form weakly bound diatomic states, crucial for the ex- 
perimental realization of the BEC-BCS crossover physics. 
In fact, it was shown in Ref. [H that the use of recombi- 
nation as an efficient way to produce weakly bound di- 
atomic molecules can be extended to systems other than 
fermionic gases. For colliding BECs at precisely-tuned 
relative velocity, the formation of molecules via 3-body 
recombination can also be used to form molecules effi- 
ciently owing to a double Bose enhancement Q . Recom- 
bination processes normally release a high amount of ki- 
netic energy, producing atomic losses that often limit the 
lifetimes of Bose- Einstein condensates (BEC) More- 
over, three-body recombination has been recognized as 
one of the most important scattering observablcs in which 
features related to the universal Efimov physics can be 
manifested Near a two-body Feshbach resonance, 

i.e., when the s-wave scattering length a is much larger 
than the range tq of the interactions, Efimov states can 
occur, causing interference and resonant effects in recom- 
bination. The experimental observation of these features 
in recombination has been recently used as evidence of 
Efimov physics in ultracold quantum gases 

From the theoretical viewpoint, quantitative calcula- 
tions of recombination for the typical alkali atoms used 
in experiments in ultracold gases are limited by the large 
number of diatomic states existing in such systems. Most 
of the available calculations for recombination for realis- 
tic systems have been confined to model systems possess- 
ing just a few bound states and/or systems with small 



scattering lengths, and even these arc challenging cal- 
culations [13. As a result, recombination calculations 
relevant for ultracold gases can only be made in the uni- 
versal regime \a\ 3> ro, by using simple potential models 
(with a few-bound states) or else by simply modeling the 
decay into all deeply bound molecular states through a 
single inelastic parameter Q. However, in ultracold gases 
experiments the condition of universality is typically not 
deeply satisfied, making it desirable to perform more real- 
istic calculations involving more sophisticated two-body 
models with, eventually, a larger number of deeply bound 
states. This paper develops a methodology within the 
hyperspherical adiabatic representation that permits the 
treatment of systems with many bound states. 

The present study still utilizes two-body potentials 
models that are, however, designed to support many 
bound states, and therefore mimic three-body collisions 
for more realistic scenarios. In the hyperspherical rep- 
resentation, the existence of many bound states leads to 
a complex set of sharp nonadiabatic avoided crossings 
in the hyperspherical potential curves at short distances. 
The large number of sharp avoided crossings creates nu- 
merical difficulties for the traditional adiabatic represen- 
tation as formulated with d/dR couplings [l^. To over- 
come these numerical difficulties, one solution is to use 
the slow variable discretiazation (SVD) method proposed 
by Tolstikhin et al [l^. The SVD method has been suc- 
cessfully applied to three-body bound-state calculations 
14 , [l^ and three-body collisions for the H -I- Ne2 system 
17|. These calculations, however, did not require numer- 
ical solution of the hyperradial equation [see Eq. ((T0| be- 
low] out to large distances. To study ultracold collision 
processes such as recombination in the large scattering 
length limit, it is crucial to solve the hyperradial equa- 
tion out to very large distances. Since application of SVD 
over the entire space is demanding in terms of memory 
and cpu-time, it is in fact much more efficient to separate 
the domain of hyperradii into two regimes. At short dis- 
tances, where many avoided crossings appear, the SVD 
method is applied, while at large distances, where the adi- 
abatic potential curves are smooth, the traditional adia- 
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batic method [T!] is utilized. 

This two-pronged strategy enables efficient calculation 
of the three-body recombination rate at low collision en- 
ergy, with extremely stable results for a two-body po- 
tential model supporting up to about 10 bound states. 
This numerical capability of calculating recombination 
with many bound states permits us to study the final 
state distribution of the recombination rate, if 3. One 
unexpectedly simple finding is that the branching ratio 
of recombination into a particular final (/) channel, de- 
fined as 

/ 

is the same for different initial (i) three-body collision 
channels. In the above equation, K^'^^^'^ is the partial 
recombination from the initial three-body channel i to a 
particular final channel /. The threshold laws for the par- 
tial recombination rates have also been considered, i.e., 
recombination processes occurring from excited three- 
body continua. These partial rates are observed to de- 
viate from the usual Wigner threshold law. Specifically, 
the energy dependence of the partial recombination rates 
display a much weaker suppression than the usual Wigner 
analysis [l^, [2^ for excited continuum channels. These 
numerical results can be interpreted as the manifesta- 
tion of a strong coupling between three-body continuum 
channels. This is further quantified through a perturba- 
tion series expansion of the scattering matrix that reveals 
the three-body recombination pathways at low collision 
energies. 

This article is organized as follows. Section II discusses 
the numerical methods adopted in this study. Section III 



shows the numerical results for the three-body recom- 
bination rates , and Section IV presents our analysis of 
the recombination pathways. Section V summarizes and 
concludes. 



II. METHOD 

The system studied here consists of three identical 
bosons with masses mi = TO2 = TO3 = m. After separat- 
ing out the center-of-mass motion, this triatomic system 
is described using modified Smith- Whitten hyperspheri- 
cal coordinates [ij, [2l| : {R, il} = {R, 9, a, /3, 7}. R is 
the hyperradius describing the overall size of the system, 
and the hyperangles and tp describe the internal mo- 
tion of the three-body system, a, /? and 7 are the usual 
Euler angles. In these hyperspherical coordinates, the 
interparticlc distances [l3| are defined as, 

ri2 = 3"l/*i?[l + sin6lsin((^-7r/6)]^/^ (2a) 

raa = S^^/^i?, [1 + sin6lsin ((^ - 57^/6)]^/^ (2b) 

rai = 3-l/^i?[l + sin6lsin(v3 + 7r/2)]^/^ (2c) 
where the i and j indices in 7'y- label particles i and j . The 
present study uses a two-body potential model potential 
in form of 

V [r,,) = Dsech" (^^^ , (3) 

where the coefficient D is negative, and its magnitude is 
chosen to be large enough to support 8 to 10 two-body 
bound states (4 to 5 s-wave bound states). 

In hyperspherical coordinates, the three-body 
Schrodinger equation is given (in atomic units, a.u.) by 



1 32 (A2 + 15/4) 



2/Z3b 2/X3bi?2 



V{R,e,^) 



il^^, (i?, ^) = E^)^, (i?, 17) , 



(4) 



where -01^' = R^^'^'^v' is the rcscalcd total wave func- 
tion. The index v' distinguishes different linearly inde- 
pendent solutions that are degenerate in energy and in- 
cludes other "exact" quantum numbers. In Eq. 
is the "grand angular momentum operator" jl3| . /^3f, = 
\J TfiYm^rfi-}, I [ra\ + m2 + mj^ ~ mj^/Z is the three-body 
reduced mass, and the total potential energy is defined 
as the pairwise sum of two-body interactions 



Y {R, 9, if) = V (ri2) + V (r23) + v (r^i) . 



(5) 



Eq. (|4]) is solved in the hyperspherical adiabatic rep- 
resentation. Like the usual adiabatic approximation, 
the hyperspherical adiabatic potentials and channel func- 
tions are defined as solutions of the following adiabatic 



eigenvalue problem: 

i/ad [R, n) (i?; n) = (R) (R; n) 



(6) 



where the adiabatic Hamiltonian, containing all angular 
dependence and interactions, is defined as 



A2 



15 



2^l3bR^ 8Ai36i?2 



V{R,9,^) 



(7) 



Therefore, the adiabatic potentials and nonadiabatic cou- 
plings, obtained by solving Eq. ([6]) for fixed values of R, 
contain all the correlations relevant to this problem. Fig- 
ures 1 and 2 show the typical adiabatic potential curves 
for the parameter m = 7.2963 x 10^, D = -5.500 x lO''^ 
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a.u. and tq = 15 a.u. Figure 1 exhibits several sharp 
nonadiabatic avoided crossings at short distances. Al- 
though in our representation sharp crossings are associ- 
ated with vanishingly small transition probabilities, these 
avoided crossings can cause several numerical difBculties 
when solving for the hyperradial motion in the traditional 
adiabatic approach. As mentioned above, such difficul- 
ties are overcome by implementing the SVD method de- 
scribed in Ref. [15}. Figure 2 shows, however, that the 
adiabatic potentials at long distances are smooth and, 
therefore, are more suitable for traditional approaches. 

The goal of our scattering study is to determine, from 
the solutions of Eq. dH, the scattering matrix S_. In 
order to accomplish such goal, the R-matrix is first com- 
puted; this ^ is a more fundamental quantity that can 
be subsequently used to determine the scattering matrix 
[see Eqs. (|35p and (jSH) below]. As usual, the i?-matrix 
TZ{R) is defined as 



'KiR) = F{R) F{R) 



(8) 



where matrices F_ and F_ are given in terms of the solu- 
tions of Eqs. dH) and © by: 



F,,, (R) - J dr!$, {fi; R)* xjj,, (ft, R) , 



(9a) 



F,,, (i?) = J dm, (17; R) (n, r) . (9b) 

In the traditional adiabatic method, ]Z{R) is calcu- 
lated by assuming an adiabatic expansion for the total 
wave function. In this case the p' — th independent wave 
function at the specified energy E is expanded in terms 
of the channel functions <I>y(i?; ft) with coefficients given 
by F,i,i{R). After projection onto the channel functions 

{R; n) from Eq. (|4]) , the resulting coupled system of or- 
dinary differential equations can be solved for this v'—th 
independent solution: 



1 



2Ai3fc dR^ 

— E 



(R) - E 



F,,, (R) 



(10) 



(R) 



F^,, (R) = 0, 



with appropriated scattering boundary conditions for the 
hyperadial solutions F_{R) and determine F_{R) in terms 
of the derivatives of F(i?) and ^,{R;fl), see Appendix 
A. The nonadiabatic couplings in Eq. ([TU]) are defined as, 

dn^AR;^)* g^^AR-^^)^ (n) 



(R) 



$^(i?;17), (12) 



are the terms that control the inelastic transitions as well 
as the width of the resonances supported by the adia- 
batic potentials Uy{R). Therefore, the accuracy of the 



results for three-body observables depends crucially on 
the accuracy of the calculated nonadiabatic couplings. 
Usually, these matrix elements are numerically evaluated 
by a simple differencing scheme, i.e.. 



$^ [R + AR: n) - $^ (i? - AR; n) 
2AR 



(13) 

This scheme works well, however, only when P,^ (R) and 
Qi^fj. (R) are smooth function of R, e.g., at large distances 
and/or for systems with just a few bound states. Clearly, 
this scheme suffers from tremendous numerical difficul- 
ties arising from sharp nonadiabatic avoided crossings at 
short distances in systems with many bound states. In 
that case, the SVD approach offers a much more stable 
and accurate approach for solving Eq. (jj). 

One key ingredient for implementing the SVD ap- 
proach is the use of the discrete variable representation 
(DVR) [H, ill. Our DVR basis functions tt; (i?) are 
defined by the Gauss-Lobatto quadrature points Xi and 
weights Wi [131. This quadrature approximates integrals 
of a function g (x) as 



g {x)dx 



N 

E 

i=i 



g {xi)w^. 



(14) 



After scaling the quadrature points and weights, the 
above equation is generalized to treat definite integrals 
over an arbitrary interval i? G [ai, 02]: 



/ giR)dR^Y.9iR^)^^ 



(15) 



where 



02-01 02 + 01 ,02-01 
m = 7^ Wi,Ri^ H . (16) 



Eq. is exact for polynomials whose degree is less 

than or equal to 27V — 1. We construct the DVR basis 
functions as 



R — Rj 
Ri — Ri ' 



which have the important property that 



TTi (Rj) 



(17) 



(18) 



Hence, matrix elements of any function H (R) obey, 

ra,2 

/ n^{R)H{R)TTj{R)dR^H{R^)S^j, (19) 

J ai 

which is usually called the DVR approximation. 

Next, the R-matrix propagation method is combined 
with the SVD approach, following the logic of Ref. [25| 
and using the DVR basis given by Eq. p7|) . For a given 
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R-matrix [Eq. ([8])] at i? = oi one uses the R-matrix prop- 
agation method to calculate the corresponding R-matrix 
at another point i? = a2, as follows. The solution -0^/ is 
expanded in the radial DVR basis tTj (i?) and in hyper- 
angles in terms of the adiabatic hyperspherical channel 
functions as 



E 



Cj^,,,n,{R)<S>,,{n;R,), (20) 



where ^^{il;Rj) is the z^— th hyperspherical adiabatic 
channel function calculated at i? = Rj. 

Substituting Eq. P0| into Eq. © yields the values of 
matrix elements of Fi,^> and F,^i,i at i? = oi and R ~ a2 
boundaries in terms of the coefficients of Eq. (|20p : 



^ vv 


(ai) = ^Cj,,yT:j (ai), 

j 


(21a) 


^ vv 


("2) = ^Cj^.v-^j (02), 

j 


(21b) 


Fvv 


(ai) ^c.j^,yOl^^TT] (ai), 


(21c) 




(21d) 


where O]^^ are 


the overlap matrix elements, 






= j dn^, in;R,)*<^>^ in;Ri). 


(22) 



Note that the determination of F_ and F_ according to the 
above expressions only depend on derivatives of the well- 
behaved DVR basis [ttj{R)] . Therefore, this approach 
is much better suited to handle the complex structure of 
avoided crossings present in systems with multiple bound 
states (see Fig. 1). 

Over an interval R € [ai, 02], the DVR approximation 
based on quadrature gives: 



TT, (R) i/ad (i?, n) TT, (R) dR « i/ad iR^, n) % . (23) 



Expansion of the Schrodinger equation in the same nu- 
merical basis functions as in Eq. (j20p and integration by 
parts yields the equation for the expansion coefficients 
Cjfj..v' (in vector notation, c^i): 



H - E 



or, equivalently. 



H -E 



(24) 



(25) 



Here the matrix elements of H and L are given by. 



Hi, 



^iVjjfi 



1 



2/i 



36 



dn, (R) diTj (R) 



dR 

1 



dR 



dR 



dR 



UfJ. 



(26) 
(27) 



Diagonalizing H over the range [01,02] gives. 



and the completeness relation of a;„. 



(28) 



(29) 



where 1 is an identity matrix. Eq. ([25]) is then rewritten 



H - e] ^ XnX^Lc,,' = ^ 



-Lc,. 



E 

n n 

(30) 

Substitution of the matrix elements of L from Eq. ([?7)) 
and insertion of the definition of F^^r and F,yi,i at ai and 
02 in Eq. (PT|) . finally gives 

(n) / \ (ri) / N. 

^ 2/i36(en-i?) 

(n) I \ (ti) I \ 

(OlX Ol g , s ,01 ^ 

~ o,. E^-^M!-' (ai),(31a) 



^li/' (02) 



E 



2Ai36 (en - £■) 

(n) I \ (n) f \ 
2/i36 (En - E) 



Ff^v' (02) 



(n) / \ (n) r \ 

Y. \}"t"' i"^^ F,.,(oO,(31b) 



nfj. 



2^36 (e„ - £■) 



where. 



(J?) = ^a;jv,„7rj (i?). 



(32) 



and Xji^^n are elements of the vector a?„. 

Our next step introduces the following matrices 



i-Kii)v, 

iRl2)v^ 

(^21).^ 
(^22).^ 



E 

71 fl. 

E 
E 
E 



In) / \ (n 
tii. (Oi) 




2Ai3h (£n — 


E) 


(n) / N (n 
Ui. (Oi) 


(«2) 


2Ai3fc (£n — 




(n) / \ (n 

ui. (02) 


(ai) 






(n) / \ (n 

ui, (02) 


(«2) 


2Ai3fc (£n - 





(33a) 
(33b) 
(33c) 
(33d) 



and after some manipulation, the matrix equation is fi- 
nally obtained that determines the R-matrix propagation 
from Ol to 02: 

^(02) =^22 -^21 B.ii+'K{ai)V^'Ki2- (34) 

In the SVD method, the overlap matrix O^*^ requires 
us to calculate the channel functions $1, (il; Rj) at every 



5 



grid point which can be very memory demanding if 
one needs to perform calculations in a broad range of 
R. At large distances, therefore, we apply the traditional 
adiabatic approach combined with the R-matrix propa- 
gation method. In the traditional adiabatic method, the 
P and Q matrixes can be calculated on a sparse grid, and 
then interpolated and/or extrapolated on a much denser 
grid and larger distances. This strategy makes the calcu- 
lation faster and it also requires less memory. The main 
difference between the traditional adiabatic approach and 
the SVD method is the use of a different three-body nu- 
merical basis. The details of this traditional approach 
and its connection with SVD method are discussed in 
Appendix A. 

Once wc have the R-matrix at large distances, the 
physical scattering matrix S_ (and its close relative, the 
reaction matrix IC) can be simply determined by applying 
asymptotic boundary conditions: 

!C^{l-rn){g-g;ny\ (35) 




-6.0x10 V ' f'-^ 1 • 1 • 1 • 1 

20 40 R 60 80 100 



FIG. 1: (Color online) Adiabatic potential curves U (R) at 
short distances 7?. Red solid lines represent the three-body 
continuum channels, i.e., the initial channels for recombina- 
tion processes, and black dashed lines represent the final re- 
combination channels. 



S = {l + ilC) {1-iicy 



(36) 



where /, /', g and g' are diagonal matrixes whose ma- 
trix elements are the energy-normalized asymptotic so- 
lutions /y, g^ and their derivatives g[, respectively. 
fi, and f/i/ are given in terms of spherical Bessel func- 
tions: UiR) = (2ii^bK/^f''' RjiAKR). gu{R) = 

1/2 

(2/i3{,fci^/7r) ' Rni^ {k^R), where and are deter- 
mined by the asymptotic behavior of the potential in 
The three-body recombination rate 

M 



2.0x10" 



Eqs. daasi) m 

is therefore given by 



K. 



1927r2(2J- 
fJ'Sbk'^ 



(37) 



where Sfi is the appropriate S-matrix elements, J is the 
total angular momentum of the system, and k ~ y'2fi3i,E 
gives the hyperradial wave numbers in the incident chan- 
nels. 



U(R) 




5000 10000 15000 20000 25000 30000 



FIG. 2: (Color online) Same as Fig. 1, but for the adiabatic 
potential curves U (R) at large distances R. This figure con- 
trasts with Fig. 1 in the characteristically smooth behavior 
at large distances. 



III. THREE-BODY RECOMBINATION RATES 

The present numerical study focuses on systems of 
three identical bosons with total angular momentum 
J = 0, with parameters adjusted to represent the "^He 
system (mi = 7712 = ma = 7.2963 x 10^ a.u. and ro = 15 
a.u.). The two-body potential depth D [see Eq. ([3])] is 
adjusted to tunc the scattering length a and explore both 
the positive- and negative-scattering length cases while 
supporting 8-10 bound states. The recombination rate 
near the unitary regime {k\a\ « 1), is explored next for 
two sets of typical parameters: D = —5.500 x 10"^ a.u. 
for the positive-scattering length case, a = 1020.36 a.u.; 
D = —5.467 X 10~^ a.u. for negative-scattering length 
case, a — —1096.07 a.u. For both cases a is much larger 



than ro (|a|/ro ~ 70) and therefore such calculations are 
solidly within the universal regime [3, [131 • 

The black dashed lines in Figs. 1 and 2 denote the 
recombination channels, i.e., the final state channels of 
the recombination process. The effective hyperradial po- 
tentials Ui, (i?) = Ui, [R) — Quv/ {2fJ.3b) for these channels 
have asymptotic behavior given by 



Uf{R) 



R^oo If {If + 1) 



E. 



(/) 

26 ' 



(38) 



where E^f, is the two-body bound state (dimer) ener- 
gies, and If is the corresponding angular momentum of 
the third particle relative to the dimer. The subscript / 
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FIG. 3: (Color online) Partial recombination rate ifg as 
a function of the collision energy E for the positive-scattering 
length case. The solid, dashed, dot, dash-dot and short 
dashed lines indicate partial recombination rates to different 
recombination channels: / — 1,2,3,4,5, respectively. The 
thick black lines, thinner red lines and thinnest blue lines in- 
dicate recombination rates from different incident channels: 
i = 1,2,3. The three green dash-dot-dot lines are propor- 
tional to ij", and E"^ as indicated in the figure. 

labels these recombination channels in ascending order, 
i.e., from high-to-low two-body bound state energies. In 
Figs. 1 and 2, red solid lines denote the three-body break- 
up channels (or entrance channels) whose asymptotic is 
described by 

~ .^oo A. (A -f 4) -f 15/4 ^ 

where A^ (A^ -I- 4) is the eigenvalue of the grand angular 
momentum operator (here. A; = 0,4,6,8..., where 
Ai = 2 is absent for symmetry considerations). The sub- 
script i labels three-body break-up channels in ascending 
order, i.e., from low-to- high eigenvalues A^. 

As Ref. [l^, [2^ points out, the asymptotic form of Ui 
determines the Wigner threshold laws for recombination, 
i.e., the low energy behavior of the recombination rate. 
A simple extension of the results of Ref. [l^ [l^] yields 
the Wigner threshold laws for all three-body channels as 

K^/^'^ (X E^^ (40) 

Our numerical results, however, show that Eq. (|40|) only 
holds for the lowest entrance channel {i = 1) while it fails 
to describe the threshold laws for higher incident chan- 
nels (i > 1). This is apparent from Figs. 3 and 5, which 
show our numerical calculations for recombination with 
positive and negative values of the scattering lengths, re- 
spectively. 

In fact. Figs. 3 and 5 illustrate that for the partial 
recombination rate from the lowest three-body incidence 
channel (Ai = 0), the threshold behavior does follow the 



FIG. 4: (Color online) Branching ratio of the calculated re- 
combination rates r^^*"^'' as functions of the collision energy 
E for the positive-scattering length case. The line-styles solid, 
dashed, dot, dash-dot and short-dashed indicate recombina- 
tion rate to different recombination channels: / = 1, 2, 3, 4, 5, 
respectively. The thick black lines, thinner red lines and 
thinnest blue lines indicate recombination rates from three 
different incident channels: i = 1,2,3. The inset shows the 
branching ratio between the deep bound states only and it 
excludes the shallowest bound state (see the text for further 
details). 



Wigner threshold law prediction, K^/^^'^ cx How- 
ever, for higher incident channels A2 = 4, A3 = 6, the 
threshold energy exponent is independent of A^ and re- 
combination rates arc only proportional to . Therefore 
the low-energy suppression for higher three-body break- 
up channels is much weaker than what Wigner 's thresh- 
old law would predict (see Eq. ([30])). Note that we have 
used different line-styles to indicate the recombination 
rate for different incident channels, and use different color 
and thickness of lines for different final channels. The 
solid, dashed, dot, dash-dot and short dashed lines indi- 
cate recombination rate to different recombination chan- 
nels: / = 1,2,3,4,5. The thick black lines, thinner red 
lines, and thinnest blue lines indicate recombination rates 
from different incident channels: i = 1,2,3. 

Another important property that has emerged from 
our numerical calculations is that the branching ratio 
[Eq. ([1])] for the three-body recombination rates into dif- 
ferent final channels are the same for the few lowest ini- 
tial channel in the low collision energy limit (see Figs. 
4 and 6). For instance, for the three different initial 
channels shown in Fig. 4, a case with positive scatter- 
ing length, the branching ratio into the highest bound 
state is about 0.35 for each of the three lowest incident 
channels throughout the energy range E < 10/j,K. Simi- 
lar results are seen for the branching ratios at negative 
scattering lengths, as is documented by Fig. 6. Note, 
however, that the branching ratios for positive scattering 
lengths are not the same for E > lO/iK (see Fig. 4) while 



10' 



10' 



FIG. 5: (Color online) Same as Fig.3 but for the negative 
scattering length case. 



FIG. 6: (Color online) Same as Fig. 4 but for the negative 
scattering length case. 



they remain the same for negative scattering length (see 
Fig. 6). This is a result of constructive interference ef- 
fects that reduces the recombination probabihty for the 
most weakly bound recombination channel for positive 
scattering length, and it is related to the universal Efi- 
mov physics 0, [201. fact, such interference effect is 
only significant in the shallowest final channel. Hence, if 
the / = 1 channel is excluded from the summation in the 
denominator on the right hand side of Eq. ([T]), the cal- 
culated branching ratio between the deep bound states 
should be the same for the whole energy range consid- 
ered, as the inset of Fig. 4 shows. 

Both the branching ratio properties uncovered in the 
present numerical exploration and the deviations from 
the recombination Wigner threshold laws can be under- 
stood using the analytical model developed in the next 
section. As we will see, these results are driven simply by 
the strong long-range coupling between the three-body 
incident channels; this analysis gives further insight into 
the pathways controlling thrcc-body recombination. 



IV. DOMINANT RECOMBINATION 
PATHWAYS 

The extensive numerical three-body recombination 
rates presented in the preceding section are now inter- 
preted in order to extract the important recombination 
pathways. Once these are identified, the surprising low- 
energy threshold behavior and the branching ratio regu- 
larities cited above become clear. 

Our model consists of carrying out a perturbation ex- 
pansion of the S-matrix and then associating each term 
to an specific pathway. As a first step, Eq. pOl) is recast 
in matrix form as 



where 



\Tr + W- El] F = 0, 



(41) 



1 



2n3b dm 



and 



IP 



dR 



QuLt 



(42) 



(43) 



The off-diagonal terms of W_ are treated perturbatively, 
suggesting that the hyperradial Green's function matrix 
should be defined as the solution of 



Tj^+H- EY\G {R, R') = 5{R-R')l 



(44) 



where U_ is the diagonal submatrix of W_, whose matrix 
elements coincide with Uy in Eqs. and ([M)) . One 

can, therefore, write the hyperradial Green's function as 



G{R,R')^-Tiif{R<)h}+^ (i?>). 



(45) 



where / and /i are both diagonal matrices. The matrix 
elements of / are the solutions of Eq. (|44| regular at 
i? = 0, and the outgoing Hankel solutions h are given by 



{R) = f{R) + ig (i?) 



(46) 



where g represent the corresponding irregular solutions. 
For the three-body break-up channels, since the cen- 
trifugal barriers are dominant, the regular and irregular 
energy-normalized basis pair fi and gi are well approxi- 
mated in terms of Bessel functions as 



9^ (R) 



\/Jj^bRJ\,+2 (kR) , 
VAW?n.+2 {kR) . 



(47a) 
(47b) 



The above hyperradial Green's function can now be used 
to expand the S-matrix in a distorted-wave Born series. 



S 



r(0) 



r(2) 



(48) 
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where 



c(0) 



Of., 



0, 



(49) 



(/ 7^ i). Ill Eq. the first order expansion of the 

scattering matrix element is simply given by, 



S 



(1) 



2TTi 



dRff (R) Wf, {R)f, (R) , (50) 



and the low energy behavior of the S-matrix elements 
can be easily determined by inspection. The integrand 
in Eq. (|50p is only significant at small values of kR where 

f, (R) = VRJx,+2 {kR) cx Therefore, 



Sf^ cx k 



Xi+2 



(51) 



In terms of the pathways, the first order S-matrix el- 
ement is the probability amplitude to transit from the 
incident channel i and then tunnel through the centrifu- 
gal barrier and scatter into recombination channels at 
short distances (i? oc tq). Therefore, if the recombina- 
tion process were solely described by this pathway, the 
low energy behavior of recombination would be given by 



K. 



1927r2(2J+l) 



S 



(52) 



recovering the usual threshold laws from Wigner's anal- 
ysis [il[23i. 

The first-order result shown in Eq. ([5^ for the low en- 
ergy behavior of recombination fails, however, to explain 
our numerical coupled-channel results [see Figs. 3 and 5] 
implying that high order perturbation terms in Eq. (|48p 
are crucial in order to determine the actual threshold 
laws. Hence we consider the second order partial-wave 
Born expansion, given by 



(53) 



where, 



POO 

h ^ dRff{R)Wf„AR)f,n(R) 

"'0 



X / dR'hl + ^{R')Wrn^{R')f^iR'), (54a) 



/•OO 

h ^ dRff{R)Wf,n{R)hi^^{R) 

X / dR'fmiR')Wm^iR')f^{R')■ (54b) 



The first integral Ii describes the quantum amplitude for 
a pathway in which the system coming inward in incident 
channel i to first scatter into an intermediate state m via 
a long-range coupling and then scatters to the final chan- 
nel / at short distances. I2 describes the amplitude for a 
different second-order pathway for which the system first 
scatters into an intermediate state m at short distances 



and then scatters into the final channel / in a second 
step. Accordingly in our analysis, the most important 
pathway for all incident channels is the one associated 
with the Ii term in Eq. (|53p. i.e., the pathways incor- 
porated in I2 are much more strongly suppressed in the 
low-energy limit. 

Interestingly, the second-order correction for the S- 
matrix element associated with the lowest three-body in- 
cidence channel {i ~ 1) can only produce a deeply-bound 
molecular channel or else an excited three-body contin- 
uum channel. In both cases, our analysis shows that 
these contributions are unimportant in the low-energy 
limit. Therefore, the threshold law for the lowest three- 
body channel is still given by Eq. ([52]) [with i = 1, 
Ai ==0]. For recombination events starting from excited 
three-body channels {i > 1), however, the situation is 
different. In this case the dominant pathway is the one 
that involves the lowest three-body continuum channel as 
an intermediate channel (m = 1), with a corresponding 
second order correction: 



^(2) 



: -2^2 / dRfi (R) Wfi (i?) /i (i?) 
Jo 

/ dR'h[+'> {R')Wu{R')f^{R')■ (55) 



The long-range coupling Wu between the lowest three- 
body break-up channel and a higher incident channel 
is dominated by the P-matrix element between the two 
channels. For different i > 1, the P-matrix element Pu 
follows the same asymptotic behavior [l^ 



Pu (R) « 



1 



(R^oo). 



(56) 



Using the above equation and definition of IVj/^ in Eq. 
(|43)) . the integral in the second line of Eq. ([55]) has the 
property that 

/•OO 

/ dR'h[+'> {R') Wu {R') h {R!) oc k. (57) 

JR. 

The integral in the first line of Eq. ([55|) is the same as 
Eq. ([5D1) . Therefore, the second-order scattering-matrix 
element for the i > 1 three-body break-up channels fol- 
lows 



f ^ (X k ^ ~^ k — k , 



(58) 



which is larger than the first-order S-matrix for channels 
i > 1 m the small k limit [see Eq. (|5ip ]. Therefore, based 
on the discussions above, the threshold behavior of the 
partial recombination rate for any incident channel can 
be summarized as 



K. 



cx 



i = l, 
i > 1, 



(59) 



which is consistent with our numerical results shown in 
Figs. 3 and 5. 
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The present analysis, therefore, demonstrates that the 
important recombination pathway for excited three-body 
incidence channels involves an intermediate transition to 
the lowest three-body incidence channel, controlled by a 
strong and long-range coupling between continuum chan- 
nels [Eq. ([56| ]. This recombination pathway also explains 
why the relative recombination rate to reach the same fi- 
nal recombination channel from different incident three- 
body channels is the same. For ultracold collisions trig- 
gered from every excited three-body incidence channel 
{i > 1), our analysis shows that the final state contribu- 
tion for recombination is mainly controlled by the cou- 
pling between the lowest three-body break-up channel 
at short distances. Therefore, the corresponding relative 
final state contribution are independent of the initial ex- 
cited three-body channel. 

V. SUMMARY 

The methodology elaborated in this paper is capable of 
calculating recombination rate and, similarly, any other 
three-body scattering observable for systems that possess 
many two-body bound states. Our numerical study was 
performed for systems with up to 10 bound states, but 
it can be extended to larger problems with a good level 
of accuracy. Although our calculations for larger systems 
might be limited by memory usage and CPU time, our 
approach still allows for the analysis of increasingly more 
complex systems. 

A key outcome is an understanding of the modified 
threshold laws for partial channel contributions to three- 
body recombination of three identical bosons with angu- 
lar momentum J = 0. Our analysis for the important re- 
combination pathways reveals that the threshold behav- 
ior of the recombination rate for excited three-body inci- 
dence channels is significantly less suppressed at low en- 
ergy than a simple generalization of the Wigner's thresh- 
old laws predicts. In addition, the branching ratio of 
recombination into any given final state / is found to be 
the same for different incident channels. 
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Appendix A: R-matrix propagation method with 
traditional adiabatic method. 

The present model described in the main text uses 
the traditional adiabatic approach combined with an R- 
matrix propagation method for large values of R, where 



the P and Q matrices are smooth functions of R. One 
advantage of using this representation is that instead of 
calculating the values of the P and Q matrices at ev- 
ery mesh point in P, we can solve the hyperangular part 
of the Hamiltonian at relatively fewer grid points. (The 
number of grid points is generally set by the characteristic 
wave length associated with the collision energy.) This, 
therefore, allows the use of interpolation and/or extrap- 
olation methods in order to generate the required much 
denser grid without memory storage problems. This Ap- 
pendix describes the approach in more detail. 
In the traditional method, V'jy' is expanded as 

V-.' (P) = ^c,,,,'^, (P) {R;n). (Al) 

A comparison of Eq. with Eq. ^ shows that 

the main differences between the traditional adiabatic 
method and the SVD method derive from using differ- 
ent three-body numerical basis sets. (Notice that in Eq. 
(|A1[) . the (P; il) are channel functions evaluated at P, 
while in Eq. the $,y {Rj;il) arc channel functions 

evaluated at Pj.) However, one can show that the expan- 
sion coefficients Cj^^^i are the same for the two expansions 
if <I>iy (P; 51) is smooth so that the DVR approximation 
can be applied, 

J dRiT, (P) (P; n) (P) « {R,; il) (A2) 

Eq. (|A2[) implies that the traditional adiabatic method 
and the SVD method are equivalent within the DVR 
approximation. Therefore, when the P and Q matrices 
change rapidly and are hard to evaluate numerically, it 
is highly advantageous to choose the SVD method; when 
the P and Q matrices are smooth, however, the tradi- 
tional method is simpler and benefits from lower memory 
storage requirements. 

Next, the details of the traditional approach are elab- 
orated and the R-matrix propagation from a point bi to 
another point 62 is explained. Insertion of Eq. (|Aip into 
Eq. (HI) yields 

J^..'(6i) = ^c,,v 7^^(61), (A3) 

j 

Pi/!.' (62) = y^^cj^^^'Wj (62), 

p..' (61) = ^[^.pP;v(6i) + ^'.m(&i)^^m^'(&i)], 

P..' (62) = 5][^.m^;./(62)+P.^(62)Pm.'(&2)]. 

As the next step, rewrite Eq. in the basis of Eq. (jAl[) 
in the matrix form of Eq. (j24p . with the matrix elements 



10 



"^fJ-Sb Jbi 

1 



tt[ (R) tt'j (R) dR5f,^ 



2/i 



36 



b2 



2A'3fc Jbi 
1 

2/X3b 



TT, (R) P^, {R) V. (R) - Tt[ (i?) P,^ {R) TTj (R) dR, 



TT, (R) S^^Wj (R) + TT, (R) P^, (R) TTj (R)] 



62 



(A4) 
(A5) 



Use of these matrix elements and replacing al(a2) with the matrix propagation. 
61(62), the same procedure as Eqs. (|28ll34p accomplishes 
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